This invention relates to a method for preparing a numerical calculation program for representing numerically physical phenomena by means of a partial differential equation describing them and displaying them schematically (hereinbelow called simulation) and in particular to a method for preparing a numerical calculation program simulating the distribution of physical quantities in a certain space such as in the analysis of electromagnetic field, heat conduction, fluid, etc.
Since such a numerical calculation program varies depending on the form of the object region for the simulation, the equation (in general partial differential equation) representing the physical law governing a physical phenomenon and conditions at the boundary of the region, it is laborious to prepare a program therefor. Consequently it is desirable to prepare automatically the program based on simple input data.
Heretofore there have been proposed some programs for such a programming (e.g. Morris, S. M. et al., "SALEM-A programming System for the Simulation of Systems Described by Partial Differential Equation", Proc. Fall Joint Computer Conf., vol. 33, pp. 353-357 (1968); Cardenas, A. F. et al. "A Language for Partial Differential Equations", Comm. ACM; vol. 13, No. 3, pp. 184-191 (1970. 3); Rice, J. R. et al. "ELLPACK progress and plans, in Elliptic Problem Solvers", pp. 135-162, Academic Press (1981); and IMSL Inc., "TWODEPEP: A finite Element Program", 3rd. ed. (1982)).
Since there are restrictions concerning the form of the partial differential equation according to the prior art techniques, which can be an object of a simulation, in order to remove them, there have been published a programming language called Differential Equation Solver and programs for preparing a numerical simulation program by using it. For example there has been published the outline of a method for preparing a program permitting numerical calculations for a partial differential equation simulating a physical phenomenon on the basis of the finite differential method in Transactions of 26th National Conference of Information Processing Society of Japan, pp. 427-432 (1983. 3).
A physical phenomenon is determined by giving a partial differential equation governing an unknown physical quantity, a spatial region which is a domain of the equation, and boundary or initial conditions therefor. A numerical calculation reproducing simulatively a physical phenomenon by using a computer can be carried out by an approximate calculation called finite differential method, when a partial differential equation describing it, boundary conditions therefor, a spatial region, instructions specifying how the spatial region is divided in a mesh form and the method how differential operators contained in the partial differential equation, the boundary conditions and the initial conditions are discretized (approximately representing the physical quantity at each mesh point governed by the partial differential equation by using the same physical quantity at mesh points surrounding it) are defined.
However, a usual numerical simulation by means of a general purpose computer can realize its function only with information stated above. Therefore, according to the prior art techniques previously stated, the control of the numerical simulation using the necessary minimum number of specifications is enabled by disposing a mechanism automatically preparing a program for controlling the general purpose computer on the basis of the above specifications.
Here a method for preparing automatically a FORTRAN program, where a partial differential equation is discretized by the finite differential method, is taken as an example. Usually, at a programming an operation is performed, by which a single module or a DO loop is formed for a plurality of blocks, for each of which the same processings can be effected.
According to the prior art techniques stated above, the program has only a function to divide the mesh points into those inside of the object region of the calculation and those on the boundary, for which a DO loop is used as a unit. For example, a rectangular object region 310 of the calculation, as indicated in FIG. 1, is divided into a plurality of regions indicated by 301 to 305 and the calculation program is prepared by forming a DO loop for each of the regions as a unit. For this reason, there was a problem that it was not possible to take measures in the case where the number of orders of the partial differentia: equation is increased, in the case where the difference rule is changed, in the case where the form of the region is not rectangular, or in the case where the partial differential equation or the difference rule is different locally. Concretely speaking, as indicated in FIG. 2, in the case where the discretization is effected with a difference rule of ##EQU1## for a second order partial differential equation, the value at a point i (point indicated by .) to be calculated is obtained by referring to the values at the points i-1, i+1 (points indicated by X) directly adjacent to the point i (FIG. 2, 311). On the other hand, for the third order partial differential equation, the value at the point is obtained by referring to the values at the points adjacent but one (FIG. 2, 312). For this reason, as indicated by 313 in FIG. 2, the reference points can be outside of the boundary, depending on the point to be calculated, and this gives rise to an inconvenience in the calculation. Further, in the case where a difference rule of ##EQU2## is used, as indicated in FIG. 3, even for a second order partial differential equation, calculation is effected while referring to the points adjacent but one to the relevant point, and this gives rise to the same problem as that for a higher order partial differential equation. When the object region of calculation is not rectangular, it is possible to form a DO loop, only if not only the mesh points are divided into those inside of the boundary and those on the boundary, but also the region, where the points inside the boundary exist, is divided, depending on the shape of the region. Consequently it was not possible to form any calculation program according to the prior art techniques.